#ifndef INTEPOLATION_H
#define INTEPOLATION_H

#include "Function.hpp"
#include <stdexcept>
#include <cmath> 
#include <vector>
#include <iostream>

class Interpolation {
protected:
    const Function& f;
public:
    Interpolation(const Function& f) : f(f) {}

    virtual double interpolate(double x) const = 0; // 纯虚函数，用于在特定点进行插值

    virtual ~Interpolation() = default;
};

class NewtonInterpolation : public Interpolation {
private:
    std::vector<double> x; // x坐标点
    std::vector<double> y; // y坐标点，存储函数值
    std::vector<double> dividedDifferences;

public:
    // 使用函数对象的构造函数
    NewtonInterpolation(const Function& f, const std::vector<double>& x)
        : Interpolation(f), x(x), y(x.size(), 0.0), dividedDifferences(y.size(), 0.0){
        for (size_t i = 0; i < x.size(); ++i) {
            y[i] = f(x[i]); // 预先计算函数值
        }

        // 初始化差商表的第一列
        for (size_t i = 0; i < y.size(); ++i) {
            dividedDifferences[i] = y[i];
        }

        // 构建差商表
        for (size_t j = 1; j < y.size(); ++j) {
            for (size_t i = y.size() - 1; i >= j; --i) {
                dividedDifferences[i] = (dividedDifferences[i] - dividedDifferences[i - 1]) / (x[i] - x[i - j]);
            }
        }

    }

    // 使用函数值数组的构造函数
    NewtonInterpolation(const std::vector<double>& x, const std::vector<double>& y)
        : Interpolation(f), x(x), y(y), dividedDifferences(y.size(), 0.0) {
        // 初始化差商表的第一列
        for (size_t i = 0; i < y.size(); ++i) {
            dividedDifferences[i] = y[i];
        }

        // 构建差商表
        for (size_t j = 1; j < y.size(); ++j) {
            for (size_t i = y.size() - 1; i >= j; --i) {
                dividedDifferences[i] = (dividedDifferences[i] - dividedDifferences[i - 1]) / (x[i] - x[i - j]);
            }
        }

        }

    double interpolate(double x_val) const {
        // 构建牛顿插值多项式
        double result = 0.0;
        for (size_t i = 0; i < dividedDifferences.size(); ++i) {
            double term = dividedDifferences[i];
            for (size_t j = 0; j < i; ++j) {
                term *= (x_val - x[j]);
            }
            result += term;
        }
        return result;
    }

};

class HermiteInterpolation : public Interpolation {
private:
    std::vector<std::vector<double>> matrix;
    std::vector<double> x;
    size_t cols;
    size_t rows;
    size_t r;
    std::vector<double> dividedDifferences;
public:
    HermiteInterpolation(const std::vector<std::vector<double>>& matrix) : Interpolation(f), matrix(matrix) {
        cols = matrix[0].size();
        rows = matrix.size();
        r = rows-1;
        x.resize(cols*r);
        for (size_t i = 0; i < cols; ++i) {
            for (size_t j = 0; j < r; ++j){
            x[i*r+j] = matrix[0][i];
            }
        }
        int fac = 1;

        dividedDifferences.resize(cols*r);
        // 初始化差商表的第一列
        for (size_t i = 0; i < cols; ++i) {
            for (size_t j = 0; j < r; ++j){
            dividedDifferences[i*r+j] = matrix[1][i];
            }
        }

        // 构建差商表
        for (size_t j = 1; j < cols*r; ++j) {
            for (size_t i = cols*r - 1; i >= j; --i) {
                if (x[i] != x[i - j]){
                    dividedDifferences[i] = (dividedDifferences[i] - dividedDifferences[i - 1]) / (x[i] - x[i - j]);
                }
                else{
                    fac *= j;
                    dividedDifferences[i] = matrix[j+1][i/r]/fac;
                }
            }
        }


    }

    double interpolate(double x_val) const {
        // 构建埃尔米特插值多项式
        double result = 0.0;
        for (size_t i = 0; i < dividedDifferences.size(); ++i) {
            double term = dividedDifferences[i];
            for (size_t j = 0; j < i; ++j) {
                term *= (x_val - x[j]);
            }
            result += term;
        }
        return result;
    }

};

class BezierInterpolation : public Interpolation {
private:
    std::vector<std::vector<double>> matrix; // 控制点矩阵
public:
    BezierInterpolation(const std::vector<std::vector<double>>& matrix) 
        : Interpolation(f), matrix(matrix) {}

    double interpolate(double x_val) const {
        return 0;
    }

    std::vector<double> interpolate_vector(double x_val) {
        size_t n = matrix.size(); // 控制点的数量
        std::vector<std::vector<double>> points(matrix); // 复制控制点

        // 逐步计算贝塞尔曲线点
        for (size_t r = 1; r <= n; ++r) {
            for (size_t i = 0; i < n - r; ++i) {
                for (size_t j = 0; j < points[i].size(); ++j) {
                    // 使用de Casteljau算法计算新的点
                    points[i][j] = (1 - x_val) * points[i][j] + x_val * points[i + 1][j];
                }
            }
        }
        return points[0]; // 返回最终计算的点
    }
};

#endif